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Abstract 



A construction of prismatic Hardy space infinite elements to discretize 
wave equations on unbounded domains S7 in Hi oc (Q), H\ oc (curl; fl) and 
-Hi oc (div; Q) is presented. As our motivation is to solve Maxwell's equa- 
tions we take care that these infinite elements fit into the discrete de 
Rham diagram, i.e. they span discrete spaces, which together with the 
exterior derivative form an exact sequence. Resonance as well as scat- 
tering problems are considered in the examples. Numerical tests indicate 
super-algebraic convergence in the number of additional unknows per de- 
gree of freedom on the coupling boundary that are required to realize the 
Dirichlet to Neumann map. 

1 Introduction 

Time-harmonic Maxwell's equations are used to model resonance and electro- 
magnetic scattering problems, that occur for example in the modelling of pho- 
tolitography processes [I] , the simulation of meta-materials [H] , photonic cavi- 
ties [55] or plasmon resonances. 

We consider the following second order elliptic equations on ft C M 3 . First 
the time-harmonic second order Maxwell system for the electric field which in 
variational form in 7Ji oc (curl; f2) is given by 



-f/ioc(curl; f2) (i? c (curl; f2)) is the space of vector fields v which are together 
with the curl V x v locally (compactly supported) in (L 2 (i7)) 3 . And second the 
Hclmholtz equation in variational form in H^ oc (il) 




Vv G H c (curl; fi). 



(1) 




(2) 
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as well as in mixed formulation for cr := ^Vm G H\ oc (div; fl) and u € Lf oc ({Y) 

/ uV ■ T - ieKcr ■ rdx = 0, W € i?(div; fi), (3a) 
Jo 

ikui; + V • crvdx = g(v) , Vu G L 2 (f2). (3b) 



n 



The constant k with positive real part is the wavenumber, g(v) contains any 
source term and e is the local permittivity or refraction index. 

We deal in this paper with two types of problems: The scattering and the 
resonance problem. The scattering or source problem consists of finding a solu- 
tion u, u or (er,u) to ([I]), ([2| or Q for a given wavenumber k > and a given 
functional g. In the resonance problem we are looking for eigenpairs to ([!]), ([2| 
or ([3]) with ,g = where n is now a complex resonance with positive real part 
and u, u or (<x,it) the resonance function. 

It is very common that ([!]) and ([2]), ^ are posed on an unbounded domain 
fi, such that radiation conditions as the Silver-Muller or Sommerfeld radiation 
condition are required to make the problem well-posed. For the numerical solu- 
tion of (JlJ or ([2| resp. ^ these radiation conditions at "infinity" are replaced 
by transparent or non-reflection boundary conditions at some artificial bound- 
ary F that separates the bounded convex computational domain f2j nt , where one 
wants to calculate a solution, from the exterior domain f2 ext [131 115) . Usually 
the exterior domain is assumed to be homogeneous without sources. For the 
presented method s is allowed to be non-constant, but there are restrictions to 



e given at the end of Section 2.3 



The ansatz to realize transparent boundary condition presented in this ar- 
ticle is based on the pole condition [331 HO] > which sloppy speaking says, that 
a solution of the Hclmholtz or timc-harmonic Maxwell equation in the exte- 
rior domain is radiating, if its Laplace transform taken along a path from the 
boundary T to infinity is holomorphic, i.e. does not have any singularities in 
the lower complex half plane. The pole condition is numerically realized by the 
Hardy-space infinite element method by mapping the lower half plane onto the 
unit disc and approximating the function there. This concept has already been 
applied successfully to the two dimensional Hclmholtz equation ([THl 1301 |3"T] ) 
and to time-dependent one dimensional Schrodinger and wave equations f|32jV 

The Hardy space infinite element method leads to tensor products of special 
functions in generalized radial direction and surface boundary functions. The 
name infinite element method refers to the classical infinite clement methods [8 , 
[TU] , where the choice of radial functions is based on the asymptotic behavior of 
Hankel functions of the first kind. There are several other approaches to realize 
transparent boundary conditions based on separable coordinates and special 
functions [M], perfectly matched layer constructions [38l El |6j, boundary 
integral approaches [21] and local high order approximations [13] . 

Most of these method have the drawback to depend non-linearly on the 
wavenumber k. Hence, for resonance problems they would lead to a non-linear 
eigenvalue problem. Although it is possible to solve the resulting non-linear 



Hardy space method 



3 



eigenvalue problems, see e.g. [25] , it is desirable to avoid them. Therefore com- 
plex scaling methods are currently the standard method for solving resonance 
problems, see e.g. [El [25]. Unfortunately these methods give rise to spurious 
resonance modes and several parameters have to be optimized for each problem. 
The presented Hardy space infinite element method has the advantage to de- 
pend linearly on k 2 and there is only one parameter and the number of degrees 
of freedom in generalized radial direction to choose. 

Apart from the treatment of the exterior problem, we use the high order finite 
elements of Schoberl and Zaglmayr [36] [39] for the bounded interior problem. 
Recent overviews of iJ(curl; Q) and 7J(div; f2) conforming method are 

[2"81 [5] and in the context of differential forms [TS] • 

The paper is organized as follows: Section[2]gives the variational formulation 
for the scalar problem ([2| as well as for ([!]) and ^ using Hardy space functions. 
In Section [3] the local basis functions for 1 (O) , i?(curl;£l), H(div, ft) and 
L 2 (fl) conforming methods are presented. In Section H the local basis functions 
are collected together and in Section [5] numerical results show the exponential 
convergence of the method. 

2 The exterior problem 

Let P be a convex polyhedron containing K \ SI. The unbounded domain f2 is 
split into a bounded interior domain f2j nt := $1 fl P and an unbounded exterior 
domain Sl cxt := O \ P. Without loss of generality we may assume that the 
surface T — Q- mt (1 Sl cx t is triangulated. In the following we first summarize 
the the basic ideas of the Hardy space method for exterior domains Jl cx t in 
one dimension and introduce our notation. Afterwards, the extension to three 
dimensions and vector-valued functions is given. 

2.1 Exterior variational formulation in one dimension 

The starting point of our method is the characterization of outgoing waves by 
the so-called pole condition [34] : A general solution to the homogeneous one 
dimensional Helmholtz equation 

- u"{x) - k 2 u(x) =0, x > 0, (4) 

is given by u(x) = C\ exp(i«;a;)-|-C2 exp(— inx), whereas u is outgoing, iff C2 = 0. 
In the original form the pole condition states that u is outgoing, if the Laplace 
transformed function U := Cu with U(s) — -Q- + -§^- has no poles with 
negative imaginary part. Hence, U has to be analytic in the lower half plane 
and belongs to the Hardy space H~(M.). 

For the definition of Hardy spaces we refer to [TT] or [19l Sec. 2.1]. Roughly 
speaking they are L 2 -boundary values of functions, which are holomorphic in 
some region. Here, we mainly use the Hardy space H~ (K) of the lower complex 
half plane {s € C | 3(s) < 0} and the Hardy space H + (S 1 ) of the complex unit 
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disk {z £ C | \z\ < 1}. A parameter dependent Mobius transform 

z + 1 



z i y s(z) — iKQ- 



1 



with no € C and 5f(rto) > is used to construct a mapping M KQ : H (M) — >■ 
which is unitary up to a constant factor by 

M Ko U(z):=U(s(z))-^. (5) 

The variant of the pole condition we are using in this paper is the following: 
u is outgoing, if U := M Ko Cue H + (S 1 ). Our aim is to derive a transformed 
variational formulation of the exterior problem, which incorporates the radiation 
condition by the choice of the Hardy space H + (S 1 ) for the transformed solution 



(see eq. (11) below). First, we need the following identity, which follows from 
basic complex analysis (see [El Lemma A.l]): For suitable functions u and v 
and their Laplace-Mobius transforms U := M. Ko Cu, V := M Ko Cv in H + (S 1 ) 
we have 



I 

Jo 



«(0«(0 = <U,V) , (6) 

/o 

where 



a(U,V):=^ [ U{z)V{z)\dz\. (7) 



The second ingredient required for the derivation of the transformed variational 
formulation is the transformation of the derivative operator in propagation 
direction. For this end we introduce the operators T ± : C x H+ (S 1 ) ^ H+ (S 1 ) 
by 

T±(uo, U)(z) := \ (u + (z± l)U(zf) . (8) 

It is easy to check, that T± are injective: If U(z) — Yl < fLo a 3 z '' ^ H + (S 1 ) and 
if the image T±{uq, U) vanishes, then |mo| = |ao| = |ai| = • ■ • , and since (oy)j 
is square summable it must be identically zero as well as uq. For solutions u 
to the homogeneous Helmholtz equation it is shown in (SUl Lemma 5.6], that 
indeed M Ko Cu belongs to the image T_(C x H + (S 1 )). 

If A4 KQ Cu = j^T-{uq,U), the transformation of spatial derivative <9j is 
given by 

M K0 C{d ( u) = T+(u 0l U) = ikoT+TZ 1 M K0 C{u) , (9) 

so we set 

df.^iKoT+TZ 1 . (10) 

This is not the only reason for using the operator T~: TZ 1 decomposes the 
function M Ko Cu into the boundary value u = u(0) and a function U with 
no contribution to the boundary value of u, which is important to ensure to 
continuity of the finite element basis functions over the boundary T (see Sec. 



3.11. 
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Figure 1: Transformation F of a right prism K to a pyramidal frustum K 



Using these formulas and definitions, the weak formulation of the one-dimensional 
Helmholtz equation ^ with e.g. Q — [— l,oo), f2 e xt = [0, oo), Qi n t = [—1,0] 
and e(x) = 1 for x € fl cx t can be transformed to 

u'v' - e k 2 uvdx + a{d^U,d^j - n 2 a(U,V) = g(v) . (11) 



For the implementation we refer to Section 3.1 or for more details to [191 Sec. 
2.4]. 



2.2 Segmentation of the exterior domain 

In three dimensions we assume that the exterior domain f2 ex t can be discretized 
by infinite non-degenerate pyramidal frustums in such a way that on each frus- 
tum e is constant and that there is one distance variable £ that allows to globally 
parameterize fl ex t by blowing up T. Such a parameterization cannot be con- 
structed locally. We refer to [3H [23l [24j |40] for details on constructing such 
a segmentation in a general two dimensional setting. These assumptions on e 
allow to treat for example layered media, infinite cones or waveguides. 

To keep the presentation simple we consider a special exterior discretization: 
Suppose that we have constructed a tetrahedral mesh of f2i n t, which induces a 
triangulation T on T. Assuming now that e is constant in fl ex t we can chose 
an arbitrary reference point Vq in P and discretize f2 ex t by infinite pyramidal 
frustums K with a triangular base T £ T and infinite faces formed by rays 
starting from Vq and proceeding through the vertices of the surface triangulation 
T: 

K := {x + £(x - Vo) £R 3 | £>0, x£T}. (12) 

£ is called a (global) generalized radial variable. To compute integrals over a 
pyramidal frustum K the parameterization 



F K {^x) :=£ + £(£- Vb) 



(13) 
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is used to transform the frustum onto a right prism K := [0, oo) x T over 
the surface triangle T, cf. Fig. [I] Using the surface gradient Vj of the two 
dimensional surface applied to the three components of x G M 3 , the Jacobian of 
the bilinear mapping Fx is given by 

J(£, x) = (x - V , (1 + Vp) g R 3x3 . (14) 

For later use it is useful to factorize 

. A o o \ 

J(^,x) = J{x) 1 + £ (15a) 

\0 o i + z) 

into on a; dependent part J(x) := (x — xq j Vp) and a £ dependent part and 
note that 

J- T (Z,x) = J(zT T | u — 11 I • d- r 'l') 




\J(£,x)\ = (l + 2 |J(i)|. (15c) 



Remark 1. In f!9f the exterior domain is the complement of a ball of radius 
a around the origin. Functions u in the exterior domain are written in polar 
coordinates with an additional scaling factor for symmetry reasons (cf. \1!A eq. 
(3.2)1): 

u cxt (r,x) = (r + lY d - i y 2 u((r + l)x). 



If the reference point Vo in (12 I is set to the origin, if P would be a ball and if we 
neglect the scaling factor, then the discretization of O ex t in this paper coincides 
with the one in 119/ . 

The canonical transformations that are used to transform basis functions v 
and their derivatives from a reference element K to basis functions v defined on 
a local element K = F{K) are summarized in the following Lemma: 

Lemma 2. Let K C M 3 such that K = F(K) with Jacobian Matrix J = Jp 
andV^ := (% v£ 1} , Vf ) T . 

1. For v G ^(K) let v o F := v. Then (Vv) oF = J~ T V^. 

2. Fori? G H(curl,K) let voF := J~ T v. Then (V X v) oF = rjr J V S ,£ X v. 

3. For v G #(div, K) let v o F := X J v. Then (V-v)of = rjr V S , S • v. 
Froo/. See e.g. [HI Sec. 3.9]. □ 
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2.3 Exterior variational formulation in if 1 (f2 e 



To extend (11) to scalar functions in three dimensions, we need to transform 



integrals over the pyramidal frustums K to integrals over the right prism K. 
Using the definitions and notations of the last section we have 



uv dx 



K 



{uoF) (voF) \J\ d£ dx. 



T JO 



(16) 



The Laplace and Mobius transform is applied in ^-direction to the functions u 
and v, i.e. we consider 



U(;x) := X KQ £( U oF(.,i))eIi-+(S 1 ), 
V(;x) := M K0 C{voF{.,x)) tH+iS 1 ) 



(17a) 
(17b) 



for x E T. In (2Ql [19] it was proven that this generalization of the one di- 
mensional pole condition leads to a radiation condition equivalent to other 
formulations of the physically relevant radiation condition. Using (|6| the in- 
finite integral J^{...)d^ is transformed into a bilinear form in the Hardy space 
H+iS 1 ): 

I uvdx= I a(VU{»,x),VV{»,x))\J(x)\dx, (18) 

where the operator T> in (jTsJ) , which is implicitly defined by 

V(M Ko Cu) =M K0 C((» + l)u) , (19) 

occurs due to the fact that the determinant of the Jacobi matrix includes the 
factor (£ + l) 2 - Straightforward calculations yield 

D 2 



(VU)(z) 



2ino 



(U)'(z) 



1 



1 



2inn 



(U) (z), 



(20) 



so T> has the following matrix representation with respect to the monomial basis 



{-■ 



,z 2 ,...}of H+iS 1 ): 



( -1 
1 



id- 



2in 



(21) 



Due to the transformation of the gradients in Lemma [2] we have 



V« • Vu dx 



K 



T JO 



J- T V^ & u o F ■ J~ T V^ & v o F \ J\d£dx. 



(22) 



We define the x dependent matrix G = (gjk)^ k=i 
to calculate 

(\J\J- 1 J- T )^,3 t ): 



|J|J-V- J and use (15 1 



5ll (x)(l+0 2 9u(x)(l + 9is(x)(l + 
g 21 (x)(l + g 22 (x) g 23 (x) \ . (23) 
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Thus (22) becomes 



Vu • Vv dx = Ag t 



® idj U 
V (id®V.g) U 



(id<gV £ )V 



(24) 



with 



.4 



G,T 



((U 1 ,...,U r ) T ,(V 1 ,...,V r f):= f 9jk{x)a{Uj{ 



,x),V k (;x))dx. (25) 



Since U £ H + (S 1 ) <g) L 2 (T), we use a tensor product notation to combine the 
operators T>, acting on H + (S 1 ) and Vj acting on L 2 (T) to operators on 
H + (S 1 ) (g) L 2 (T) (see e.g. [29]). id denotes the identity operator. Note, that 



the matrix G in (25) can be replaced with |J| in order to get a short notation 
for the right hand side of ( 18 ) as well. 

Similar to ( 11 ) a variational formulation of ([2| using the Hardy space H + (S 1 ) 
consist of the integral over the bounded interior domain Q- ln t and the integral 
over the unbounded exterior domain Q cx t, which is given by the sum over all 



surface triangles of the integrals (18) and (24): 



Vu • Vv — en uvdx 



TGT 



A, 



G,T 



(id®V £ )£7 




2?<% ® id J V 



A 



e\J\,T 



(U,V) 



(26) 

The coefficient function e in can be incorporated into the exterior integrals, if it 
depends for each pyramidal frustum only on the surface variable x. Numerically, 
the integrals over the surface triangle T will be approximated using a quadrature 
formula, while the bilinear-form a can be computed analytically as presented in 



2.4 Exterior variational formulation in if (curl; f2 ext ) 

Since a solution u to ([lj satisfies the vector valued Helmholtz equation and 
since the Silver-Muller radiation condition is equivalent to the Sommerfeld ra- 
diation condition for the Cartesian components of u (see [7]), we can extend 
the pole condition to exterior Maxwell problems using it for each component of 
u = (tii, U2, U3) T ■ Hence, with the transformation of _ff(curl;£!) functions of 
Lemma [2] we define 

/ M KQ Cm o F(;x) \ 
U(.,£) J T (x) \vM K0 £u 2 oF(;x) , x e T. (27) 
\VM K0 Cu 3 oF(;x)J 

The transformation of the mass integral in iJ(curl, K) has already been used 



for gradient fields and is given in analogy to (24) by 



IK 



u ■ vdx = Aqt ({(V®id)U u U 2 , U 3 ) T , {{V ® id)V u V 2 , V 3 f) . (28) 
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Due to the transformation of the curl operator in Lemma [2] we define C 

tek)lk=i :-|ir 1 i T iwith 



1 



J T J{^x) 



I c u (i) ci 2 (x) c 13 (x) 

I (1+C) 2 l+£ l+« 

°fff c 22 (x) c 23 (x) 

V c 32 (x) c 33 (x) 



The factor (£ + 1) 1 in (29) leads to the integral operator I 
Hardy space formulation of the curl curl integral: 

VxU'Vxvdi = 



(29) 

V^ 1 in the 



A" 



A 



C,T 



(x®vf)u 3 -{x®vf)u 2 \ 

' id(g)V (2)Ni TT - f^^-^\TT. 



Ui - 9f ® id )U a 



V (d e ® id J J7 2 - (idsVj^J Ui J 



( 



V 



■ ®vi 1} ) V3- (i<8>vi 2) )y 2 \ 



id®V 



% ® id V 2 - id«)VV 



7(1) 





\ 


) v 3 




)v lJ 


) 



(30) 

As in the previous section, the integrals in ( 28 ) and ( 30 ) form the exterior part 
of@: 



ff(v) 



V x u • V x v — en 2 u ■ vdx + ( 

T<£T 



(((x^)u 3 -(x^)u 2 \ 



A 



C,T 



id(g)V? } ) Ui - [d 6 Old) J7 3 



V 



d £ ®id t7 : 



idOV^ ) J7i 



/ 



V 



id®V 



d £ <g> id V2 



7(1) 



)v 2 \ 


\ 


) v 3 




)Vi) 


) 



(31) 



- A S Q t x (((2?®id)^i, C/ 2 , *7 3 ) T , ((23® id)V^, V 2 , ^) T ) 

2.5 Exterior variational formulation in if (div; Q ext ) 

The construction follows analogously to the H (curl; ft) functions of the previous 
section using the transformation of i7(div; Q) function of Lemma [2] 



'VVM Ko C UioF(»,x)' 
U(.,z) := \J(x)\J^(x) ( VM K0 C u 2 oF(;x) 

VM K0 C u 3 oF(;x) 



x£T, 



(32) 



for u = (ui,u 2 ,u 3 ) T £ H(div,K). Hence, the mass integral in i7(div, K) is the 
integral in (30 1. The stiffness integral follows by straightforward calculations 
using the transformation of the divergence operator. 
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3 Local exact sequence of tensor product spaces 

For a bounded domain C R 3 , which is diffeomorphic to the unit ball the 
sequence 

C -A H\n) A #(curl,fi) ^> £T(div,n) A L 2 {Q) A {0} (33) 

is exact, i.e. jC = ker(V) (where j maps a E C to the constant function Q — > C, 
x i-> a), Vif^fi) = ker(Vx), V x H(curl,ft) = ker(V-) and V • F(div,Q) = 



L 2 {Q). On the two dimensional surface Y the sequence ( [33] ) decomposes into 
two exact sequences 

C H X {Q) ^ i?(curl,fi) ^ L 2 (ft) A {0}, (34a) 

C A H 1 ^) % JT(div,fi) ^ L 2 (n) A {0}, (34b) 

with the rotated surface gradient operator V^r := (— V~ 2 \ V^) T and the scalar 
curl operator curl^ = V^r • . In the following • _L always denotes the rotation of 
a 2d vector, vector space or operator with values in C 2 . 

For Maxwell eigenvalue problems, due to exact sequence property there exists 
an eigenvalue with an infinite dimensional eigenspace consisting of all gradients 
of H 1 functions. If the finite clement method is not constructed carefully, the 
numerical results will be polluted by artificial eigenvalues coming from this 
eigenspace (see e.g. [31 Sec. 6.2]). Hence, following [551 Chapt. 5] or [35] we 
try to carry over the exact sequence property of the continuous spaces to the 
discrete spaces. Until now we cannot prove a discrete de Rham diagram with 
commuting projection and differential operators for the discrete spaces we are 
going to construct, but at least they build an exact sequence locally. 

In the interior domain we use a discretization with tetrahedral elements and 
the high order elements of [39, Chapt. 5.2.6] with the local exact sequence 

C -A P p A (P'P- 1 ) 3 ^ (PP- 2 ) 3 A (PP- 3 ) 3 A {0}. (35) 

In the following we focus on the elements for the exterior domain and the cou- 
pling on the artificial boundary. 

There, the pyramidal frustums K are mapped into right prisms K (see 
Fig. [lj. Therefore the local elements are build by tensor products of one- 
dimensional infinite elements for £ 6 [0, oo) and triangular elements on the 
2d surface element T consisting of subspaces Wt C 7? 1 (T), Vt C H (curl, T) 
and Xt C L 2 (T). The construction of the local sequence is motivated by a 
modified tensor product of two chain complexes [171 Sec. 5.7]: Combining the 



two surface sequences ( 34 ) we use the cochain complex array 



W ( <g> W T — W ( ® V T ^ ® X T 

<9 5 ®idj 9 ? ®id- L | 8 t ®id| (36) 

W'^Wt W' t ®V£ W'^Xt. 
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The tensor product chain complex is given by the sums over the diagonals: 



W% (g) W T 



(i 



w 4 ( 



o 

■idigiVi 



id(g>V^ 
^(giid 1 



Xi 



( 9 4 <giid id®V 4 - ) 



(37) 



WL ®X T . 



gence operator. Hence, (37 1 discretizes in some sense 



If r C K , the operators in (37) are the standard gradient, curl and diver- 



A H (curl, K) — 



H(div,K) L 2 {K). 



(38) 



The term ,,in some sense" is used, because the finite element basis function 
will not belong to iJ 1 (iC), H(curl,K), H(div,K) and L 2 (K), since the radial 
parts of them belong to the Hardy space H + (S 1 ). The definition of the correct 
function spaces would be very complicated. In [HO Sec. 3.2] the function space 
is given for scalar functions with spherical artificial boundary. Here, we con- 
fine ourselves to the discrete spaces, but we indicate in the following with the 
notation 

V:=(%V«vi 2) ) r , (39) 
that d(: is an operator of a subset of H + (S 1 ) into H + (S 1 ). 



3.1 iJ 1 -element 

We briefly discuss the iJ 1 -elements proposed in [191131]. For each segment K 
we use the local tensor product space 

W k = W 6 <g> W T , (40) 

whereas Wt is the surface element of the surface triangle T of the H 1 tetrahedral 
volume element and W(_ a space in radial direction. Wt equals the usual poly- 
nomial space P P (T) for the two-dimensional triangular element and therefore it 
has the dimension |(p + 2){p + 1). 

In radial direction we have to discretize the Hardy space H + (S 1 ), where a 
L 2 -orthogonal basis is given by monomials. Using the first N + 1 monomials 
IIjv := span{z°, ...,z N } and the operator 7~_ of ([8]), we define the discrete space 

:= T-(C x Hn) with the basis functions 

*-i := — r_(l,0) and^- := — T-(0, j = 0,...,N. (41) 

Since (£ _1 M£ %) (0) = for j = 0, ...,N and (C 1 M^ (0) = 1, the 
basis function is used to couple the Hardy space infinite elements of the 
exterior domain O ext to the finite elements of the interior domain Q- mt . 

Let us collect all basis functions for the prism K and assign them to a 
vertex Vi (i — 1,2,3) of the surface triangle T, an edge Eij between vertex 
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^---Q ~-4>*,- (giu) 1 ^ 



Figure 2: basis functions in Wj> = W% ® Wr 

Vi and Vj , the surface triangle T, an infinite ray Ri corresponding to a vertex 
Vi, an infinite face corresponding to the edge or the prism K itself. 
Following [SHI Chapt. 5.2.3] a basis function w of Wt can be a 

1. vertex basis function w Vi , 

2. edge basis function wf' J , I = 1, ...,p — 1, or a 



3. element basis function 



1 = 1,... i(p-2)(p-l). 



Hence, the basis functions of W^- = <8> Wt sketched in Fig. [2] are 
1. vertex functions W Vi := <g> u^, 



2. edge functions W, 



*_i <g> w, w , / = 1, ...,p - 1, 



3. surface functions Wf := \I>_i ® w ; T , 2 = 1, |(p — 2)(p — 1), 



4. ray functions W, 



Hi 



to**, fc = 0,...,iV, 



5. infinite face functions W, 



Ft) ._ 
k.l 



*k®iuf M , fc = 0, ...,JV, Z = 1, ...,p — 1 and 



6. segment functions := $ k ®wf, k = 0, ...,iV, Z = 1, |(p-2)(p- 1). 

The degrees of freedom for the segment basis functions are interior degrees of 
freedom, while the others have to coincide between neighboring infinite segments 
and on the surface with the tetrahedron in f2i n t. In total, there are 



1 



dimW^ = -(p + 2)(p + l)(J\r + 2) 



(42) 



degrees of freedom. 



3.2 i/(curl)-element 



Due to the sequence (37) we use for the approximation of elements in H{c\iv\) 
the space 

W>®W T \ 

W i ® V T J ' { ' 



V k = 



Hardy space method 



13 




Figure 3: basis functions of V% 



with the surface finite element space Vp of the iJ(curl) volume element of order 
p — 1. As in the previous case Vr is the vectorial triangular element space with 
two components, i.e. Vr = (P p_1 ) with dimVr = {p+ l)p- For (10) leads 
to the basis functions 

rl>-x := T+(l, 0) and ^ := T+(0, (.)''), j = 0, ...,N. (44) 

If v ; B * J € Vt denotes the edge basis functions of the i?(curl) surface trian- 
gular element and vf the surface triangle basis functions, the basis functions in 
Vfr are 

1. edge functions vf' 3 := ( ^ E ] , I = 1, ...,p, 

\W_i ® Vj / 

2. surface functions VT := | T ^ T ] , I = 1, (p — 2)p, 



3. ray functions Vf* := ^ rK J, jfe = -1, ...,iV, 

4. two types of infinite face functions 



Vi 



(a) V^ 1 : " v , J,* V./ 1 p, 

(b) V^ 2 := ^ ^ ),* = -!, iV, i = 1, ...,p - 1 
5. and two types of segment functions 

( a ) V W : = i v r ) , k = 0, N, I = 1, (p - 2)p and 







'fc » vf 




(b) Vff := (>* ® «f ') . fc = _ x JVT, / = 1, .... |(p - 2)(p - 1). 
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See Fig. [3] for a scheme of the curl-conforming basis functions. Note, that 
only the tangential directions indicated by the arrows are continuous over the 
boundaries. The segment basis functions have vanishing tangential components 
on the edges, rays, faces and the surface and are therefore interior basis func- 
tions. The tangential components of the other functions vanish on the boundary 
parts to which they are not assigned. For example the tangential component 
ipk ® w Vl of a ray function on R\ vanishes due to the vertex function w Vl on 
the rays R2 and R3 as well as on the infinite face F23. It does not vanish on 
the surface, but there it is the normal component, which does not need to be 
continuous for H(curV) functions. The tangential component on the surface is 
the second component of the vector and therefore zero on the whole segment. 

If we link together the degrees of freedom for neighboring elements, we get 
a i?(curl)-conforming method with 

dimV^ = (JV + 2)Q(p + 2)(p + l) + (p+l) P y (45) 
degrees of freedom for each segment. 

Remark 1. It is not essential to use for the basis functions ipk- The special 
choice of the basis functions was necessary to divide into boundary degrees 
of freedom and interior ones, which is useful in order to guarantee continuity 
of the global basis functions at the segment boundaries and hence to get a con- 
forming method. But as we have seen for the ray function on R\, there is no 
coupling between the ray and face basis functions and the basis functions in the 
tetrahedron. 

Moreover, all of the functions £ _1 M.'^ ipk have non-vanishing boundary val- 
ues. Since T+(C,Hn) t = IIv+i we could also directly take {z°, z N+1 } as 
basis functions. The only reason for us using ipk is, that in this way the basis 
functions in Wi are exactly the derivatives of the basis functions in . 



3.3 #(div)-element 

Starting from V^, we calculate 





id®Vjr 



id®V^- 
d e ® id- 1 



(46) 



The space H(div$) is due to pil a rotated if(curl £ ) and V%W T C V?. More- 
over, the normal components of the basis functions on each face of a iZ(div) 
tetrahedral element include the scalar curl^-fields of Vt- Therefore we chose 



Qk 



K 



W 6 ®Q T 



(47) 



where Qt = P' p ~ 2 is the finite element space for the normal component on the 
surface of a iJ(div) volume element of order (p — 2) and takes the place of Xt 
in the third term of (37 1. 
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Let q T be a basis functions in Qt- Then, the basis functions of consist 

of 

1. surface functions Qf := ^ Ql ^ ) / = ^ \p{p — 1), 

F ( \ 

2. face functions Q,*, 3 := / e- A 1 - , k = — 1, iV, / = 1, p 

3. and two types of segment functions 

(a) Q^ 1 := ( * fe ® «H , fc = 0, N, I = 1, ±p(p - 1) and 



(b) : =U*®(v?') ± J' ,fe = " 1, -' JV ' Z = 1 '-' (P " 2)P - 

The basis functions of ® Qt ensure the continuity of the normal component 
on the boundary surface and those of <g> V^- the continuity of the normal 
component on the infinite faces. 
For each segment there are 







dimQ^ = (iV + 2) Qp(p_l) + (p+l)p^ 



(48) 



degrees of freedom. 



3.4 L 2 -element 

The L 2 -element has no coupling on the boundaries. Nevertheless, we compute 
the divergence of Q^: 

*•«*-( ) ■ ( $ I % ) = * ® * + ^ ® (49) 

Both parts fit together and therefore we define 

X k ^W' i ®X T , (50) 

with a two-dimensional triangular element Xt — P p ~ 2 of order p — 2. In the 
interior domain we use volume elements of order p — 3, but since L 2 -elements 
do not require continuity, there is no conflict. 

dimX k = ±p(p-l)(N + 2). (51) 



4 Properties of the sequence 

The analysis for this sequence is far from being complete. However, we are able 
to prove the exact sequence property locally. 
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4.1 Local properties 

By construction it holds for each infinite element 

VW k C {ve^|Vxv = 0} ) 
VxV k C {v G | V • v = 0}, 
V-Q k C X k . 

Note, that constant functions are not outgoing. Hence idC = ker(V) is not 
meaningful in this context. By counting the degrees of freedom we get the 
following theorem. 

Theorem 1. The discrete tensor product spaces defined above build a local exact 
sequence, i.e. 

W K A V K ^ Qk A X K A {0}. (52) 

Proof. For a linear operator T on a finite dimensional space D we have the 
identity dimD = dimT(D) + dimker(T). From (49) we conclude, that dim V • 
{Qk) — \v{v — + 2), since d r W r <g> Q t has this dimension. With (51) we 
get V • (Q k ) = X k and with Q dimker(V-) = (N + 2)(p + l)p. Again, since 
dim V x V k > (N + 2)(p + l)p, we deduce V x V K = ker(V-) and using ( |45| 
dimker(Vx) = |(p+2)(p+l)(iV+2) . Last, dimVW k > |(p+2)(p+l)(iV+2) 
and hence VW k = ker(Vx) and ker(V) = {0}. □ 

4.2 Global properties 

For the global finite element spaces in the exterior domain we define the spaces 
for the different kinds of degrees of freedom 

Wj xt := |J {W v }, (53a) 
v 

:= |J {Wf | I = l,...,p-l}, (53b) 

E 

W£t := IJ W T I f = l,-,|CP-2)Cp-l)}, (53c) 

T 

: = U i W k \k = 0,...,N}, (53d) 
i?, 

:= IJ {<! I fc = 0,...,JV,l = l,...,p-l}, (53e) 

W^cxt := IJ {^W I * = 0,...,iV,i = l > .. .,i(p_2)(p-l) (53f) 

K 

and finally collect them together to 

W cxt := U U U W* t U < xt U W* t . (54) 
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Figure 4: Hclmholtz resonator 

The spaces V ex t, Qext and X ext are denned analogously. In this way theorem [I] 
should carry over into the global finite clement spaces. 

5 Numerical results 

Although there is no convergence result of the Hardy space method for elec- 
tromagnetic scattering or resonance problems, we can expect from the results 
for scalar problems (see [T5] and the convergence plots in [3T]) super algebraic 
convergence of the method w.r.t. the number of degrees of freedom in the Hardy 
space. In order to substantiate this predication, we give in the following some 
numerical tests for scattering as well as resonance problems. Especially the re- 
sults of the challenging problems in Sections |5.2[ |5.4| and |5.5| indicate, that only 
a few basis functions in the Hardy space are needed. 

For all computations the mesh generator netgen [35] together with the high 
order finite element code ngsolve is used for the interior problem. Using a source 
term g(v) =^ in or ^ and a given wavenumber k > leads to a system 
of linear equations, which we solve directly with the package PARDISO [33]. A 
resonance problem consist of finding k with 5f(re) > and non-trivial resonance 
functions u solving the homogeneous (no sources, vanishing boundary condi- 
tions) problems 0, ^ or ([3| and results in a generalized matrix eigenvalue 
problem of the form 

Sun = K 2 h Muh 

with symmetric, non-hermitian complex matrices S and M. This problem is 
solved using a shifted Arnoldi algorithm. 
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Figure 5: Real part of three resonance functions and corresponding resonances 
of a single Helmholtz resonator with volume 11, length of the neck 5cm and cross 
sectional area of the neck 3cm 2 . 

5.1 Resonances of a Helmholtz resonator 

In the first example we compute resonances of a three dimensional acoustic 
Helmholtz resonator and show the system response to a source given by an 
incident plane wave Ui{x) = e lKd ' x with d = (— 1/v2j 0, — 1/V2) T and fixed 
wave number k > 0. A single resonator consists of a ball with radius 0.062, 
which is connected via a cylinder with radius 0.0098 and length 0.5 to the half 
space {i£l 3 : x 3 > 0}, cf. Fig. [Zj The infinite half space is truncated using 
a half ball in the case of the scattering problem and a cube in the case of the 
resonance problem. 

We assume sound-hard boundary conditions at the boundary of the resonator 
as well as at the infinite plane {x £ K 3 : x 3 — 0}. For the Hardy space method 
no sources outside the finite element domain fii nt are allowed and the solution 
of the problem needs to be outgoing in f2 0Xt . Therefore, in fii nt we compute the 
total field u and in £7 cxt the field u s = u — Ui—u r with u r being the reflected wave 
of the unperturbed half space u r {x) = e lKd ' x with d = (— 1/^/2, 0, 1/V2) T . u s is 
outgoing with vanishing Neumann boundary values at the infinite plane, since 
it is the perturbation of Ui + u r caused by the resonator. Using u in Q int and u s 
in Q ext leads to a jump condition at the artificial boundary and an additional 
boundary term for the jump in the Neumann values. 

Fig. [5] gives cross-sections of a single Helmholtz resonator with three differ- 
ent resonance functions and the corresponding resonances. For c = 343m/s we 
calculate the frequency / = c5R(k)/27t 115.8Hz, of the most relevant first res- 
onance which has a quality factor Q = 3?(k)/|9(k)| 631. The approximative 
formula (see e.g. [37]) for a Helmholtz resonator with volume V, neck length L 
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Figure 6: Real part of the total field to an incident plane wave Ui{x) = e" 
with d = (-1/V2,0,-1/V2) T and k = 2.1241. The values of the total field 
vary between ±2000, but the figure is limited to values between ±10. 



and cross sectional area of the neck S gives / = g^ry j^y ~ 133Hz. 

For these computations a mesh of 2477 tetrahedrons with finite element order 
7, the number of degrees of freedom in radial direction N = 15 and ko = 8 lead 
to approximately 155000 degrees of freedom. Note, that the third resonance in 
Fig. [5] has multiplicity 2. 

Fig. [6] shows the real part of the total field to an incident plane wave with 
wavenumber k = 2.1241 for 9 uniformly distributed resonators. For this system 
of resonators the first 9 resonance frequencies are between 114.8Hz (k sa 2.103) 
and 116.1Hz (k m 2.126) with quality factors between from 162 to 864765. For 
the resonance computations of the system of resonators we used a mesh with 
24394 tetrahedrons, finite element order 7, N — 15 and Kq = 4. 

5.2 Convergence test of a magnetic dipole 

In this example we resolve a magnetic dipole located at a point y 



E y (x) = V x x 



e iK\x—y\ 

\x - y\ 




in two different geometries (see Fig. [7]). E y is a radiating solution of ([lj with 
given wavenumber n > and an unbounded domain which does not contain 
the center y. 

In the upper part of Fig. [7] the interior domain is the intersection of two 
balls with radius 5 and 6 ( Q- ia t — B 6 \ B 5 ), k — 1 and the center of the dipole is 
the origin (y — 0). For the interior boundary dB 5 we use a Dirichlct boundary 
condition given by the tangential part of E y . For the exterior boundary DBq 
we have three different cases: First we again use the exact Dirichlct boundary 
condition in order to compute the error of the finite element discretization of 
r2i n t with polynomial order 6. Second we use the first order absorbing boundary 
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Figure 7: Right: H (curl, f2i nt )-error of the HSM w.r.t. the number N of degrees 
of freedom in radial direction compared to the finite element error and the error 
of a first order absorbing boundary condition for the two different domains Q- ln t 
on the left; Left: Cross-section of one Cartesian component of a magnetic dipole 
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condition VxEx v—inE — coming from the Silver-Muller radiation condition, 
which can be applied without spending additional degrees of freedom. Last we 
use the Hardy space method with ko = 10, reference point Vb = (0, 0, 0) and 
N = 0, . . . , 8 , which leads to 102279 (N = 0) up to 229975 (N = 8) degrees of 
freedom. 

In the lower part of Fig. [7]the interior domain is the intersection of two cubes 
ft int = [-3.2, 3.2] 3 \ [-3, 3] 3 , k = 5 and the dipole is located at y = (1, 1.5, -1) T . 
The finite element method needs 555480 degrees of freedom for polynomial order 
5; together with the Hardy space method with ko = 12.5 and Vb = (0, 0, 0) T we 
get 584642 up to 1080374 degrees of freedom. 

Both cases show a fast convergence of the Hardy space method, such that 
setting TV = 4 (N = 6) suffices to reach the finite element error. 

5.3 Cavity resonances 

Here we search for resonances k ^ and radiating electric fields E ^ solving 
for ft = M 3 \ K and K = [-1.2, 1.2] 3 \ ([-1, l] 3 U [1,1.2] x [-0.2, 0.2] 2 ) . Ad- 
ditionally, E has to satisfy at dK the perfectly conducting boundary conditions 
E x v — with v being the outward normal vector. This problem is an extension 
of the two dimensional acoustic open cavity problem in [19 . A similar acoustic 
problem is treated in |27j with boundary element methods. 

In Fig. [8] the absolute value of two resonance functions on a cross-section 
of the interior domain ftj nt = [— 1.7, 1.7] 3 n ft is shown. For a closed cavity 
(ft = [—1, l] 3 ), the resonances are positive and analytically given by (see [T]) 
K = sfl + m + for l,m,n£ No such that Im + In + mn > The resonance 
functions in Fig. [8] belong to resonances close to the second cavity resonance 
k = V3f with multiplicity 2. 

Fig. [9] shows the real and imaginary part of the computed resonances for 
two different discretizations. For the first we use the domain ft; nt = B2.5 (1 ft 
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Figure 9: Resonances of an open cavity for two different discretizations: First 
^int = -B2.5 H fl, kq — 5, N — 10 and finite element order 5. Second Qi nt = 
[—1.7, 1.7] 3 fl fi, Kq = 4 — i, N = 6 and finite element order 6 



with kq = 5, N — 10 and finite element order 5 and in total 358924 degrees 
of freedom. For the second dicretization fi int = [— 1.7, 1.7] 3 fl fl, n Q = 4 — i, 
N = 6 and finite element order 6 lead to 390548 degrees of freedom. Both 
discretizations give similar results for the cavity resonances near the real axis. 
The multiplicity of the resonances is not visible in Fig. [9] ft is the same as 
expected from the resonances of the closed cavity given in pQ. The exterior 
resonances with in absolute values larger imaginary parts are mostly identical 
for the two discretizations, but for the resonances at the bottom of Fig. [9] the 
discretizations are too coarse. 



5.4 Resonances of GaAs pyramidal micro-cavities 

A second example of cavity resonances is taken from |22) . The cavity is a 
pyramid with height 0.14142 and a quadratic base of length 0.28284 which is 
turned up-side down and mounted on top of an infinite pyramid. Choosing the 
apex of the pyramids as reference point Vq — (0,0,0), the infinite pyramid is 
bounded by the infinite rays in direction (1, 1, — 1) T , (—1, 1, — 1) T , (—1, 1, — 1) T 
and (— 1,— 1,— 1) T . The computational domain f2j nt is a cuboid given by the 
vertices (-0.2,-0.2,-0.1) and (0.2,0.2,0.2). The pyramids are made of GaAs 
(e = 12.25). In contrast to the first example the exterior domain outside the 
plotted interior domain consist of two different materials, namely air (e = 1) 
and the infinite GaAs pyramid (e = 12.25). 



In Fig. 10 a cross-section of the field intensity of the resonance field for the 
resonance k « 12.3034 — 0.8816i is given. 

To compute a reference solution f2j nt is discretized by 2477 tetrahedrons 
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Figure 10: Left: Cross-section of the field intensity of the resonance function 
to the resonance k ss 12.3034 — 0.8816i. Right: Relative error of the computed 
resonance w.r.t. the number of degrees of freedom in radial direction for different 
finite element orders. 
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Figure 11: computed and exact acoustic resonances with O = R 3 \ B\, Vt mt = 
[—1.2, 1.2] 3 \ Bi, k q — 5 — i, reference point Vq — and N = 15 



and finite elements of order 7 leading to approximately 1.5 million unknowns 
for the interior domain. For the exterior domain 18 degrees of freedom in the 
radial direction cause in total 2 million unknowns. The error in Fig. [lO] is the 
relative error of the computed resonance with respect to the reference resonance 
k w 12.3034 — 0.8816i for different finite element orders and different numbers 
of degrees of freedom in radial direction. The results indicate that depending 
on the finite element order only 8 to 10 degrees of freedom in radial direction 
are needed. 



5.5 (H(div), L 2 ) resonance test 

In the last numerical test we compute acoustic resonances k outside a sphere 
with radius 1. In this case the exact resonances are given by the roots of 
the spherical Hankel functions of the first kind, see [501 Example 3.24]. The 
multiplicity of a resonance is In + 1 if the resonance is the root of the nth 
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Figure 12: Left: Cross-section of one of the resonance functions to the resonance 
k ~ 1.754 — 1.839i; right: relative error for different finite element orders and 
different iV of the computed resonance against the analytical one with f2; n t = 
[—1.2, 1.2] 3 \ Bt, kq = 7 — 0.3i and reference point Vq = 



Hankel function. Instead of solving the Hclmholtz problem ([2]), we solve the 



mixed formulation ^ and use the iJ(div) and L 2 elements of Sections 3.3 
and [331 

In Fig. 11 the computed resonances n for Q,- lnt = [— 1.2,1. 2] 3 \ B\, finite 



element order 3, k$ = 5 — i, reference point Vq = and N = 15 are compared 
to the analytical ones. Again, for the resonances with in absolute values larger 
imaginary part the discretization with in total 397316 degrees of freedom is too 
coarse. 



Fig. 12 shows the relative error of one computed resonance against one of 



the roots of the third spherical Hankel function of the first kind. In comparison 
to the results of Sec. |5.2| and Sec. |5.4| more degrees of freedom in radial direction 
are needed, since the Hardy space method has to resolve the Hankel function. 
In total 27256 (order 1, N = 0) up to 549920 (order 4, N = 14) degrees of 
freedom are used. 
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